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Abstract 

A symplectic, symmetric, second-order scheme is constructed for particle evolution in 
a time- dependent field with a fixed spatial step. The scheme is implemented in one 
space dimension and tested, showing excellent adequacy to experiment analysis. 
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1. Introduction 

Particle motion in a space-time dependent field is a classical fundamental problem of 
dynamics. It is generally well solved in most settings under most kinds of requirements. 
However, most solutions deal with formulations of the dynamics where the independent 
variable is time, viz. they aim at computing a function y such that the motion reads 
x = y{t). Yet in some settings one actually describes motion by the reciprocal function 
t = y' 1 , so that the motion reads t = t(x). 

One such instance is the propagation of electrons in a traveling wave tube, where 
it is natural to record particles when they pass at a fixed probe location, instead 
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of getting a snapshot of their locations at a given time. Similar physical contexts 
are met in other particle beam devices, such as accelerators, klystrons, free electron 
lasers, electronic tubes for wave amplification.. To some extent, this description 
is somewhat analogous to eulerian descriptions of flows in hydrodynamics. 

If one were interested in the evolution of a single particle, one could merely compute 
its motion as x = y(t) and deduce its "schedule" function t = r(x) and related quan- 
tities as functions of spatial position. In this respect, many symplectic methods are 



available (see e.g. Refs [12|, |15| for an overview), especially for separable hamiltonians 
of the form H(p, x, t) = K(p) + V(x, t). However, to describe a beam of many particles 
during their spatial progression, it is reasonable to follow them consistently in space, 
to generate numerical data sampled at the same (possibly many) space positions. It 
then becomes awkward to first evolve them in time and afterwards reconstruct their 
progression in spatial terms by interpolations. 

For this purpose we reformulate in Section [2] the particle equations of motion, using 
the streaming variable x as independent variable (see figure [1]). Since the original 
particle dynamics is hamiltonian, we ensure that the new description be symplectic 
by first expressing the action principle in terms of the timetable function r. In the 
corresponding hamiltonian picture, the variable conjugate to r is the energy (, and the 
generator of motion is momentum V. 

In Section [3] we stress our requirements on the scheme and consider alternative 
strategies. Then we construct a first order symplectic scheme for the particle motion. 
The implicit part of the step can be performed either through algebraic solution of a 
cubic equation, or through a Newton iteration : we compare both procedures. Next 
we construct the adjoint, first order symplectic scheme, which also requires a Newton 
iteration, and we check its accuracy. Finally, we combine the direct and adjoint schemes 
to obtain a second order symmetric, symplectic, fixed Ax scheme. 

In Section H] we benchmark our algorithm by analysing the particle motion in the 
field of a single harmonic wave, viz. we solve the pendulum motion in a galilean frame. 
Numerical simulations for realistic beam data generate beam deformations shown in 
Section 

Section O focuses on the evolution of a beam launched in presence of two harmonic 
waves. Simulations are confronted with experimental observations of the beam col- 
lected at the device outlet. Special attention is paid to the reproduction of a devil 
staircase structure, characteristic of the chaotic behaviour of the system, taking into 
account the finiteness of the experimental device. 

In summary, experimental data often relate to limited interaction times, while 



2 



numerical evidences and theoretical discussions of chaotic dynamics often deal with 
trajectories followed for long times in a compact domain of phase space. The agreement 
of our simulations with experimental evidence assesses the relevance of our algorithm 
to such experimental settings. 

2. Evolution with respect to space 

Rewriting the equations of motion in hamiltonian form with respect to space is 
straightforward in the symplectic formalism (see Section [2T2T) . However one may wish 
first a more pedestrian derivation, from the classical action principle. 

2.1. Lagrangian viewpoint 

The action for a non-relativistic particle with mass m moving along a one-dimensional 
axis Ox in a time dependent potential V(x,t) reads 



where y is a continuously different iable function of time t G [to,ti] C R, subject to the 
constraints y(t ) = x and y{t\) = x\, and the dot denotes derivative with respect to 
t. The lagrangian is 



In the following we restrict all trajectories to the class of strictly monotone, increasing 
functions, viz. inf 4o<t<tl y(t) > 0. For these functions the reciprocal function r : 
[xo,Xi] — > [to,ti] : x i — y t(x) exists ; r is unique and also strictly monotone, increasing, 
continuously differentiable with 




to 



(1) 




(2) 




(3) 



where the prime denotes derivative with respect to x. 

To rewrite dTJ as a space-integral, we introduce the new lagrangian 




(4) 



u 



so that 



5[r;xo,xi,to,^i 



]:=S[t 1 ;t ,ti,x ,x 1 ] = / C(r(x),T'(x),x)dx . (5) 
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It is convenient to introduce the opposite to the (usual definition of) canonical mo- 
mentum conjugate to r, 

and to perform the Legendre transform of —C, defining 

V{(,t,x) :=Cr' + C(t,r',x) (7) 

so that in the new variables the canonical Hamilton equations read 

dr dV 

dx d( 



dC dV 
dx dr 



(8) 
(9) 



For the classical lagrangian ([2]) the new variables are the usual energy and linear 
momentum, 

C = ^js + V(x,t(x)), (10) 
V = y/2m(C-V(x,T(x))). (11) 

2.2. Hamiltonian viewpoint 

The hamiltonian formulation of dynamics provides a direct path to the latter equa- 
tions. Indeed it suffices to consider the symplectic 2-form 

dtu := dpdx - dHdt (12) 

where p = my is conjugate to x and H = ( is conjugate to t — r, and to consider x as 
the independent variable along trajectories instead of t. The minus sign we introduced 
in ©-([Tj) ensures the usual signs in du. 

Of course, if the potential does not depend explicitly on time, ( is a first integral. 

The above requirement, that y nowhere vanishes, will be strengthened below by 
imposing V > p m i Q > for some p m i n to ensure appropriate numerical accuracy. 

3. Discrete time canonical transformations 

To preserve the symplectic structure while integrating (|TU|) - ([TT|) we use a sequence 
of canonical transformations. Note that (fTTT) does not allow using a splitting method 



such as leap-frog [12j, [15|], as it is not the sum of integrable generators. 
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Figure 1: Sketch of the algorithm principle in (x,t) plot. The continued line is a trajectory (identical 
for both schemes, with small steps). Crosses sample positions (every 50 At) as obtained by leap-frog 
with constant time step, circles sample times (every 50A:c) as obtained by our approach with constant 
space step. Straight lines are guides to the eye. 



In order to sample the spatial evolution regularly (as desired e.g. to follow many 
particles in parallel), we use a fixed spatial step Ax, a strategy recommended e.g. by 
Henon [l3| to obtain simple and accurate Poincare maps (for a space periodic potential, 
the step Ax is best taken as a fraction of the wavelength). A fixed spatial step also 
avoids generating noisy-like "spurious frequencies" in the simulated dynamics. Note 
that we could also use a time-integrator (leap-frog or other) with an adapting time 
step, At = mAx/V [13j], but it is not manifestly symplectic as V is not constant. 
Another approach, using a fixed time step and interpolating the trajectories for the 
spatial mesh, would raise the issue of constructing a symplectic interpolation scheme 
(of desired order). In this work we settle for the manifestly symplectic, fixed step 
approach, which can process many particles in parallel. 

Our first integrator '■ (t, (, x) i-> (f, (, x + Ax) is chosen to provide an explicit, 
first-order approximation for the time increment. It is generated by 



so that the system 



F(t,() = t( + V(tX,x)Ax 

C = ( + d T V(r,(,x)Ax, 
f = r + d ( V(T,(,x)Ax, 



(13) 

(14) 
(15) 



is a first-order, symplectic approximation to (JS}-©. As the second equation is explicit 
with respect to time f , the first equation is implicit with respect to the energy (. For 
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the special case of momentum (TIT]) , it actually leads to a cubic equation, 

~ ~ ~ 777 

C 3 - (V + 2C)C 2 + (2V( + C 2 )C - VC 2 - -(d T V) 2 Ax 2 = , (16) 

which can be solved algebraically for (. Here V and d T V are computed at (x,t). 
Equation ffl6|) usually has three real solutions, two of which being Ax-close to ( : the 
relevant root is such that (£ — Q d T V > 0. 

It is advantageous to express (}T1]) in the form C = C + (C — V)®, so that 

a = -a(l + ar 1/2 (17) 

with the single parameter 

a = y^/2((-V)- 3/2 d T VAx (18) 

calculated at (x,r). The fixed point equation fflTj) is easily solved by the Newton 
method, which selects the "good" cubic root automatically (as aa < 0) and converges 
very fast, especially if a is small. Analytically, this method stresses the small parameter 
a controlling the accuracy of our scheme : it involves a balance of the potential evolution 
d T V and the spatial step Ax against (( — V), viz. against the particle velocity Vjvn. 
For small velocity, the algorithm deteriorates - at worst it will miss turning points 
where V changes sign (which is forbidden by our assumptions on trajectories in the 
action principle). 

Let Z%j. : (t,C, x ) ^ and Z^ x : (r, (,x) h-> (r, £, x) denote respectively 

the cubic and Newton solvers. For perfectly accurate computations, they coincide and 
may be denoted identically Z^ x . Note that Z^ x is not symplectic, as 

det DZ Ax (t, C, x) = [1 + d c d T V(r, (, x)Ax]" 1 . (19) 

With the new energy (, ffl5l) immediately provides the new time, defining the map 
Tax : (t, (,x) i-> (t, C,x). This map is not symplectic either, as 

det DT Ax (t, C, x) = 1 + d ( d T V(r, (, x)Ax . (20) 

Finally, we advance position, with Z& x : (r,C,x) !->■ (r,(,i + Ai). The resulting 
integration scheme J-Ax = %Ax ° Ta^ ° Z& x is symplectic by construction, within ma- 
chine accuracy, as the planar map Ta x ° Z& x is area-preserving : det D{Ta x ° Z Ax ) = 
(det DTax) (det DZ Ax) = 1- It is first order only. 
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Table 1: The largest single step error for five particles with different initial velocity launched in the 
field of a single wave with A = 0.1, k = 0.2, v$ = 1, 4> = 7r /4- 

The variables advanced with (TT4"]) and (fl5|) are in the form 

t" 

f £ t + r'Ax + —Ax 2 + 0(Ax 3 ) , (21) 
C = C + C'Ax + ^-Ax 2 + 0(Ax 3 ), (22) 

where r' and (' have been approximated with —d T V(r, (, x) and d^V(r, (, x). It follows 
that the most influential theoretical error, in every step of integration, is given by 
t"Ax/2 for f and ("Ax/2 for (. Table [T] compares the upper estimated theoretical 
errors, related to r" and and the maximum real simulation errors for five particle 
initial velocities. Slower particles are found affected by larger errors as expected. 



To obtain a second order, symmetric scheme, we consider the adjoint map [12 
which is also symplectic, generated by the function 

F*(t, = fC - V(f, (,x + Ax) Ax (23) 

so that 

T = f - d c V(f, C, x + Ax) Ax , (24) 

C = (-d T V(f,(,x + Ax) Ax. (25) 

For momentum fill I) , both equations involve the potential V(x + Ax, f), which implies 
that one first solves (J24J) with respect to the new time f by a Newton algorithm, and 
then computes the new energy £ by (125]) . This defines the symplectic map = 
Z lx ° Tlx ° ^Ax ■ (t, C, x) ^ (f, C, x + Ax). 

One easily checks that T^ x is self-adjoint, while Z\ x o Z_Ax and T£ x o T-Ax reduce 
to identity up to machine numerical tolerance (typically 10~ 15 ). 
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Figure 2: Discrepancy, between direct evolution Jai and backward evolution J-l Ax , for energy (left 
panel) and arrival time (centre) as functions of position, for five particles in the field of a single wave 
with A — 0.1, k — 0.2, ity = 1, (f> = tt. Particles injected at the origin at t = 0, with velocities 1.1, 
1.5, 2, 2.5 and 3.5, move as shown on right panel. 

Finally the composition 

•^Acc = ^Ax/2 ° Fax/2 = 2 Ax/2 ° VHx/2 ° ^Ax ° Tax/2 ° Zax/2 (26) 



is its own adjoint. It is thus symmetric and therefore second order 12 



4. Validation : particle dynamics in a single wave 

We test our schemes with the time dependent potential of a wave 



V(x,t) = Acos(kx-kvt + (p) (27) 

where A, k,v,cj) are respectively the amplitude, wavevector, phase velocity and phase 
of the wave. Rescaling energy (and amplitude), space and time enables one to set 
m, k and v to unity, and the choice of the origin of time or space eliminates 0. Its 
integrability makes this dynamics a good benchmark. 

The accuracy of the determination of the adjoint map is checked by iterating first 
Tax for Ax = 0.01 from x = 0tox = L = 300, and then T*_a x from x = L to x = 0, 
for five particles. Figures [2] display the discrepancies A£ and At as functions of x for 
each particle and confirm that J 7 * Ax = J 7 ^ to numerical accuracy. The order of the 
algorithms and their accuracy is further analysed in figure [3j comparing the first order 
and second order schemes for the motion of a particle with initial velocity v- m = 1.5 in 
the field of a wave with A — 0.1, = tt, k — 0.2, v $ = 1 over a length L = 100. 
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Figure 3: Numerical accuracy for the first order algorithm (red, upper lines) and second order (blue, 
lower lines) versus calculation time per spatial unit length. Darker colour for the Newton method 
solution to the cubic equation, lighter colour for polynomial solver (slightly slower than Newton). 



5. Beam dynamics in a single wave 

A second accuracy check is provided by the Poincare section of the beam by the 
positions imodi = 0. The return map for (r, () variables is symplectic. As the 
particle motion is integrable (it reduces to the pendulum by a Galileo transformation 
to the reference frame comoving with the wave), each orbit must generate section 
points on lines satisfying the algebraic relation 



C = — ± v^2mH - 2mAcos[k(x - v^t) + 0] (28) 

where H is a constant. In particular, the motion on the wave separatrix corresponds to 
H = A. For H > A, this relation defines two branches for all times r, which correspond 
to faster or slower circulating particles, while for H < A the relation defines the upper 
and lower part of trapped motion inside the wave's cat eye. Figure H] shows that 
numerical trajectories perfectly reproduce these lines. 

To assess the relevance of the algorithm to experiment we also follow the deforma- 
tion of a beam of electrons injected in a single wave, e.g. in a traveling wave tube. As 
particles are accelerated or decelerated by the wave, the beam velocity profile is de- 
formed. We thus inject N particles at x — 0, equally distributed over one time period 
of the wave, and plot the histogram of particle velocities as a function of abscissa x. 

In figure [5] a cold beam is made of particles injected resonantly with the wave 
velocity, = v$ = 1, with m = k — 1. The wave amplitude, A = 0.002, determines 



Figure 4: (left) Poincare section in (r, £) variables, of the return map after one wavelength L. Trajecto- 
ries for five trapped particles injected with v- m = 1.8, 1.85, 2, 2.1, 2.25, for five untrapped particles with 
«in = 1.4, 1.5, 1.6, 1.7, 2.3, and for one particle injected on each separatrix branch. Wave parameters 
A = 0.02, v,j, = 2, k = 0.2, <j) = 7r ; particle mass m = 1. (right) Exact section lines f|2S[) . 



the bounce frequency Wb = so that particle oscillations in the wave trough have a 
spatial period Lb = 2nv l j > /u!\> = 140.5. The particles bounce indeed and, in agreement 
with the rotating bar approximation [13], most of them reconvene every Lb/2. Only 
the ones injected at times close to (4> + 2nn) / (kv^) (for integer n) enter the wave close 
to the X point and follow closely the inner side of the cat eye separatrix. The time 
these particles need to overcome half a wavelength can be arbitrarily long, so that they 
mark the boundary of the wave resonant domain at velocities ± 2^/ A/m. 

For particles injected with a velocity outside the wave cat's eye, the beam is mod- 
ulated. A significant difference between (x, v) plots for propagating beams and the 
more familiar (x,v) Poincare sections of particle-in-wave dynamics (see e.g. 11, 10, 9[) 
is the asymmetry between faster and slower particles, obvious in Figure [6j 

In particular, for particles injected with a velocity above the wave cat's eye, the 
beam is moderately asymmetric. But for a particle injection velocity within the capture 
range [v^ — 2^J 'A/m, + 2y/A/m], the picture gets strongly deformed, as part of the 
beam is trapped as in figure [5] while part of it moves outside the wave cat's eye. 



6. Particle dynamics in two waves : resonance overlap and chaos 

The motion of a particle in the field of two waves is a paradigm of hamiltonian 
chaos. In our formulation, Poincare sections are given by xmodL = 0, and the return 
map is symplectic, hence area-preserving in conjugate variables (r, (). Figure [7] displays 
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Figure 5: Velocity distribution function of particles along the x axis, injected at the wave phase 
velocity, v\, = v^, = 1, with m = k = 1 and wave amplitude A — 0.002. 




Figure 6: Velocity distribution function of particles along the x axis, when injected at Vb = 1.6 (upper 
left), «b = 1.5 (upper right), v\> — 1.4 (lower left) and = 1.3 (lower right), above the wave phase 
velocity, = 1, with m = k = 1 and wave amplitude A = 0.04. The cat's eye boundaries lie at 
^0 + 2VJ = 1.4 and i> - 2\[A = 0.6. 
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Figure 7: Poincare section in (r, Q variables, of the return map after one wavelength L, for 15 particles, 
with mass m — 1. Wave parameters are k\ = 1, fc 2 = 1/2, w^i = 1, — 2, (fix — 4>2 — ft- Wave 
amplitudes A\ — A-i yield overlap parameter s = 0.5 (left), 0.66 (centre) and 1 (right). 



this Poincare section, showing the growth of the chaotic domain for increasing wave 



amplitudes, and the destruction of KAM tori 11 



The corresponding transition to large scale chaos by increasing the resonance over- 
lap parameter s = 2(a/Ai + / \\^<t>2 — %i|-\/m) is also observed by recording the 
particle velocities at a fixed traveled distance L , after being injected at a fixed velocity 
fi n . As seen in figures El a cold beam injected in a wave cat's eye spreads over the 
velocity interval spanned by this cat eye, and if the beam is injected outside cat eyes 
it remains confined between the velocities of KAM tori on either side. Beam velocity 
spreading (also called heating) has been used to diagnose resonance overlap, and our 
numerical scheme is compared with experimental data js[ in figure [HJ 

Moreover, the transition to large scale chaos in phase space is known to occur step- 
wise. For increasing wave amplitudes, successive KAM tori get destroyed, so that the 
beam invades velocity domains resulting from the merging of capture regions corre- 



sponding to "secondary" resonances [11(. The accessible velocity interval for the beam 
injected in one wave then grows like a devil's staircase, the higher steps corresponding 
to the merging with major secondary resonances. Figure [9] compares these domains 
obtained both numerically and experimentally 16, Q, Isf] . While experimental data 
are blurred due to recording accuracy, numerical data have limited resolution due to 
the large number of particles (here only 25000) needed for the sharp observation of a 
threshold. Nevertheless, the agreement is quantitatively satisfactory. 
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Figure 8: Velocities of particles after interaction with two waves, versus injection velocity, (top) 
Experimental data with dotted lines marking cat eyes boundaries ; (bottom) numerical results, (left) 
Non-overlapping cat eyes, s = 0.63 ; (right) overlapping trapping domains, s = 1.5. 
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Figure 9: Velocities of particles after interaction with two waves, versus amplitude of waves A\ = A%. 
Cold beam injected at phase velocity of one wave, (left) Numerical results for Vh — ityi = 1.1487, 
v<p 2 = 0.8609, ki = 0.6529, k 2 = 1.7424. (right) Experimental data. 
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7. Conclusion 



The scheme has proved its relevance to describe accurately particle motion in a 
given field, along a single space dimension. It also provides new pictures of known 
behaviours, which complement more familiar, usually more symmetric plots in (x, v) 
space. 

Extending our approach to three space dimensions is rather straightforward pro- 
vided the particles stream along a given coordinate, say x. Higher-order symplectic 
schemes can also be constructed by composing several maps T 7i Ax an d F*. Ax with ap- 



propriate substeps 7j [12|]. More challenging is the issue raised by the particle feedback 
on the wave field, calling for a self-consistent space-based model of particles and waves 
evolution, in the spirit of models used for weak plasma turbulence jsj. 

It will also be interesting to apply this scheme to model the many-waves regime of 
weak plasma turbulence, where particle velocity undergoes a chaotic transport over a 
wide range 0, Q, H, 0, 0] , as properties of this transport are still controversial. 
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Appendix A. Explicit derivation of ( 117ft 

To reduce (fT6|) to ( IT7|) . note that the derivatives of momentum ( ITT]) read 

d c V(r,(,x) = mV' 1 , (A.l) 
d T V(r,(,x) = -mP- l d T V , (A.2) 

and consider the dimensionless (both physically and numerically relevant) quantity 
a = (( — ()/(( — V), which characterizes the changes in particle energy per step and 
must be small for an accurate calculation. 

Equation ( JTTj) reads, with arguments (r, (, x) in V and (r, x) in V, 

a = -d T VAx/(C -V) (A.3) 

^drVAxiC-Vr^iC-V)- 1 (A.4) 
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1/2 



(A.5) 



(A.6) 



using the dimensionless a denned by ({TBI . The fixed point equation (|T7|) then follows. 

Appendix B. Explicit flowchart of the second order algorithm 

The symmetric scheme (126]) may be expanded as follows: 



• (r, C, x) ' y (r, C, x) : solve C = C + d T V(r, (*, x) Ax/2 for C* ; 

• (r, C, x) ^ (t*,C, x): t* = t + «9 c P(t, (*, x) Ax/2 ; 

• (r*, (*, x) y-¥ (r*, (*, : x = x + Ax ; 

• (t*, C*, 5) (t 1 , C ^) : sc, l ve r* = f - d^V(f, C*, x)Ax/2 for f ; 

• (f, C, x) (f , C, x) : C = C - d r V(f, C, x)Ax/2. 
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